{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"../common/rfsoc_book_banner.jpg\" alt=\"University of Strathclyde\" align=\"left\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Notebook Set B\n",
    "\n",
    "---\n",
    "\n",
    "## 02 - Quantisation\n",
    "In this notebook we will investigate quantisation, which is the process of mapping an analogue amplitude acquired by an ADC to a discrete set of values. We will explore the quantisation process and discuss resolution and discrete levels. We will briefly discuss the total achievable dynamic range for a given $N$-bit converter and then investigate errors and noise that may occur during quantisation.\n",
    "\n",
    "## Table of Contents\n",
    "* [1. Introduction](#introduction)\n",
    "* [2. Quantisation](#quantisation)\n",
    "    * [2.1. Resolution and Levels](#resolution-and-levels)\n",
    "    * [2.2. Dynamic Range](#dynamic-range)\n",
    "    * [2.3. Noise](#noise)\n",
    "* [3. Conclusion](#conclusion)\n",
    "\n",
    "## Revision\n",
    "* **v1.0** | 21/01/23 | *First Revision*\n",
    "\n",
    "---\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Introduction <a class=\"anchor\" id=\"introduction\"></a>\n",
    "Before we explore quantisation, let us begin by generating a test waveform. We are unable to generate a truly continuous waveform to be discretely sampled and quantised. This is because we are working on a digital system that is unable to generate continuous waves. Instead, we will use a high sampling rate to approximate a continuous wave."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import rfsoc_book.helper_functions as hf\n",
    "\n",
    "fs_cont = 10e3 #Hz\n",
    "ts_cont = 1/fs_cont\n",
    "n_cont  = np.arange(0, 2000)\n",
    "fd = 13 #Hz\n",
    "\n",
    "x_cont = n_cont/fs_cont\n",
    "y_cont = np.sin(2*np.pi*fd*ts_cont*n_cont)\n",
    "\n",
    "hf.plot_timeseries(\"Continuous Waveform\",\n",
    "                   [x_cont], [y_cont],\n",
    "                   ['continuous'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Quantisation <a class=\"anchor\" id=\"quantisation\"></a>\n",
    "We investigated the sampling process in the [previous notebook](01_sampling.ipynb), where an analogue waveform was converted to a digital waveform using a regular sampling period. Quantisation is the process of mapping each sample of the digital wave to a discrete set of possible amplitude levels.\n",
    "\n",
    "### 2.1. Resolution and Levels <a class=\"anchor\" id=\"resolution-and-levels\"></a>\n",
    "An ADC performs quantisation using $N$ bits of resolution, which means it can represent $2^N$ different possible values. A linear quantiser (shown in Figure 1) maps the samples of an input wave to a discrete set of quantisation levels.\n",
    "\n",
    "<figure>\n",
    "<img src=\"./images/linear_quantisers.png\" style=\"width: 50%; vertical-align: middle;\"/>\n",
    "    <figcaption><b>Figure 1: A plot demonstrating the functionality of a linear quantiser.</b></figcaption>\n",
    "</figure>\n",
    "\n",
    "We can obtain the quantisation step $q$ of an ADC as\n",
    "\n",
    "$$\n",
    "q = \\frac{V_{\\text{max}}}{2^{N-1}}.\n",
    "$$\n",
    "\n",
    "For a 4 bit ADC, with $V_\\text{max} = 1$, the quantisation levels are given below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Vmax = 1\n",
    "N  = 4\n",
    "q  = Vmax/(2**(N-1))\n",
    "levels = np.arange(-1, 1, q)\n",
    "\n",
    "print('Quantisation levels for 4 bits:')\n",
    "for idx, level in enumerate(levels):\n",
    "    print('Decimal: ', idx-2**(N-1), ' |  Discrete Level: ', level)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can discretely sample the continuous waveform from before without quantisation, as shown in the following code cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set sampling parameters\n",
    "fs = 100 #Hz\n",
    "ts = 1/fs\n",
    "n  = np.arange(0, 21)\n",
    "\n",
    "# Perform sampling process\n",
    "x = n/fs\n",
    "y = np.sin(2*np.pi*fd*ts*n)\n",
    "\n",
    "# And plot\n",
    "hf.plot_timeseries(\"Discrete Signal Sampled at 100Hz, No Amplitude Quantisation\",\n",
    "                   [x_cont, x], [y_cont, y],\n",
    "                   ['continuous', 'discrete'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The wave can be quantised using 4 bits of resolution, which results in the plot given below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Perform quantisation (4 bits)\n",
    "yq4 = np.zeros(np.size(y))\n",
    "N  = 4\n",
    "q4  = 1/(2**(N-1))\n",
    "levels = np.arange(-1, 1, q4)\n",
    "for idx, sample in enumerate(y):\n",
    "    # https://stackoverflow.com/questions/12141150\n",
    "    yq4[idx] = min(levels, key=lambda x:abs(x-sample))\n",
    "\n",
    "hf.plot_timeseries(\"Quantisation using 4 Bits\",\n",
    "                   [x_cont, x], [y_cont, yq4],\n",
    "                   ['continuous', 'discrete'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that many of the sample points do not lie directly on the original continuous waveform. Increasing the resolution to 6 bits improves the quantisation step (and accuracy) of the quantisation process, as shown below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Perform quantisation (6 bits)\n",
    "yq6 = np.zeros(np.size(y))\n",
    "N  = 6\n",
    "q6  = 1/(2**(N-1))\n",
    "levels = np.arange(-1, 1, q6)\n",
    "for idx, sample in enumerate(y):\n",
    "    yq6[idx] = min(levels, key=lambda x:abs(x-sample))\n",
    "\n",
    "hf.plot_timeseries(\"Quantisation using 6 Bits\",\n",
    "                   [x_cont, x], [y_cont, yq6],\n",
    "                   ['continuous', 'discrete'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2. Dynamic Range <a class=\"anchor\" id=\"dynamic-range\"></a>\n",
    "We can compute the Dynamic Range of an $N$-bit converter as follows.\n",
    "\n",
    "$$\n",
    "\\text{Dynamic Range} = 20\\text{log}_{10}(2^N)=20N\\text{log}_{10}(2)=6.02N\n",
    "$$\n",
    "\n",
    "Therefore, the dynamic range of an ADC should increase as the number of bits increase. We can demonstrate this relationship in the following code cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N_vec = np.arange(2, 17, 1)\n",
    "DR_vec = N_vec * 6.02\n",
    "for idx, DR in enumerate(DR_vec):\n",
    "    print('Number of bits: ', N_vec[idx], ' | Dynamic Range (dB): ', np.round(DR))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.3. Noise <a class=\"anchor\" id=\"noise\"></a>\n",
    "Quantisation noise can be modelled as Additive White Gaussian Noise (AWGN). All quantised samples are subject to an error. The range of possible error values are in the range $-q/2$ to $q/2$ if we assume the input signal to the quantiser is random. A Probability Density Function (PDF) denoted as $p(e)$ can be formed as given in Figure 2.\n",
    "\n",
    "<figure>\n",
    "<img src=\"./images/q_error_PDF.png\" style=\"width: 25%; vertical-align: middle;\"/>\n",
    "    <figcaption><b>Figure 2: A plot of the PDF of quantisation noise for the quantisation interval $q$.</b></figcaption>\n",
    "</figure>\n",
    "\n",
    "A technique known as dithering can be used to break-up periodic noise components in a quantised signal. Before we implement dithering, we will first generate a test waveform."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set sampling parameters\n",
    "fs = 100 #Hz\n",
    "ts = 1/fs\n",
    "L  = 2048\n",
    "n  = np.arange(0, L)\n",
    "fd = 10.001\n",
    "\n",
    "# Perform sampling process\n",
    "x = n/fs\n",
    "y = 0.7071*np.sin(2*np.pi*fd*ts*n)\n",
    "\n",
    "# And plot\n",
    "hf.plot_timeseries(\"Discrete Signal Sampled at 100Hz, No Amplitude Quantisation\",\n",
    "                   [x[0:100]], [y[0:100]],\n",
    "                   ['discrete'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The log-scale power spectra of this waveform can be inspected below. Notice that the primary lobe is over the desired sine wave frequency of 10Hz. There are no other significant tones in the sampled wave."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Y = np.fft.fft(y)\n",
    "Y_norm = np.abs(Y)*2/L\n",
    "Y_log = 20*np.log10(Y_norm)\n",
    "\n",
    "fig = plt.figure(figsize=(10,4))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(n[0:L//2+1]*fs/L, Y_log[0:L//2+1])\n",
    "axes.set_title('Power Spectra of the Sampled Wave')\n",
    "axes.set_xlabel('Frequency (Hz)')\n",
    "axes.set_ylabel('Power Spectra (dB)')\n",
    "axes.grid(True, 'Major')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will now quantise the waveform above using 4 bits of precision."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "yq4 = np.zeros(np.size(y))\n",
    "N   = 4\n",
    "q4  = 1/(2**(N-1))\n",
    "levels = np.arange(-1, 1, q4)\n",
    "for idx, sample in enumerate(y):\n",
    "    # https://stackoverflow.com/questions/12141150\n",
    "    yq4[idx] = min(levels, key=lambda x:abs(x-sample))\n",
    "    \n",
    "hf.plot_timeseries(\"Quantisation using 4 Bits\",\n",
    "                   [x[0:100]], [yq4[0:100]],\n",
    "                   ['discrete'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Errors have been introduced due to quantisation. We can easily inspect spurious tones introduced by quantisation by plotting the power spectra of the quantised wave."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "Yq4 = np.fft.fft(yq4)\n",
    "Yq4_norm = np.abs(Yq4)*2/L\n",
    "Yq4_log = 20*np.log10(Yq4_norm)\n",
    "\n",
    "fig = plt.figure(figsize=(10,4))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(n[0:L//2+1]*fs/L, Yq4_log[0:L//2+1])\n",
    "axes.set_title('Power Spectra of the Quantised Wave')\n",
    "axes.set_xlabel('Frequency (Hz)')\n",
    "axes.set_ylabel('Power Spectra (dB)')\n",
    "axes.grid(True, 'Major')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that the power spectra of the quantised wave contains tones that were not in the original sampled wave.\n",
    "\n",
    "We can derive the power of the quantisation error. The square of the error signal, weighted by the probability of error, is integrated across all possible error values. We can perform a definite integral as the range of possible error values is between $-q/2$ and $q/2$.\n",
    "\n",
    "$$\n",
    "n_{\\text{ADC}} = \\int_{-\\inf}^{\\inf}e^2p(e)de = \\int_{-q/2}^{q/2}e^2p(e)de\n",
    "$$\n",
    "\n",
    "Since $p(e)=1/q$ for all values of $e$, the error power is\n",
    "\n",
    "$$\n",
    "n_{\\text{ADC}} = \\frac{1}{3q}e^3 \\Big|_{-q/2}^{q/2} = \\frac{q^2}{12}.\n",
    "$$\n",
    "\n",
    "We can compute the noise power below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "noise_power = (q4**2)/12\n",
    "noise_power"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A normal (Gaussian) distribution can be used to generate a dither signal for the given noise power above. For a Gaussian random variable, the average power is calculated by summing the mean value and the variance, which is also known as the second moment.\n",
    "\n",
    "$$\n",
    "\\text{average power} = \\mu^2 + \\sigma^2\n",
    "$$\n",
    "\n",
    "The expected mean value of our dither signal is $\\mu = 0$. Therefore, the average power now becomes,\n",
    "\n",
    "$$\n",
    "\\text{average power} = \\sigma^2.\n",
    "$$\n",
    "\n",
    "The average power of our dither signal will be $n_{\\text{ADC}}$. Now we obtain,\n",
    "\n",
    "$$\n",
    "n_{\\text{ADC}} = \\sigma^2.\n",
    "$$\n",
    "\n",
    "We can generate a normal distribution to act as our dither signal. The Numpy random module has a normal distribution method that accepts a mean value and standard deviation $\\sigma$. We can obtain the standard deviation by simply computing the square root of the average power, as below.\n",
    "\n",
    "$$\n",
    "\\sqrt{n_{\\text{ADC}}} = \\sigma.\n",
    "$$\n",
    "\n",
    "Now, running the code cell below, we generate our dither signal using `np.random.normal(mean, standard deviation, length)`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dither = np.random.normal(0,np.sqrt(noise_power),L)\n",
    "\n",
    "hf.plot_timeseries(\"Dither Signal\",\n",
    "                   [x[0:100]], [dither[0:100]],\n",
    "                   ['discrete'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The dither signal is added to the original sampled signal before it is quantised, as below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "yd = y+dither\n",
    "\n",
    "yq4d = np.zeros(np.size(y))\n",
    "for idx, sample in enumerate(yd):\n",
    "    # https://stackoverflow.com/questions/12141150\n",
    "    yq4d[idx] = min(levels, key=lambda x:abs(x-sample))\n",
    "\n",
    "hf.plot_timeseries(\"Quantised Signal with added Noise\",\n",
    "                   [x[0:100]], [yq4d[0:100]],\n",
    "                   ['discrete'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can observe the spectrum of the quantised signal with dither and overlay it on the same plot as the original quantised signal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Yq4d = np.fft.fft(yq4d)\n",
    "Yq4d_norm = np.abs(Yq4d)*2/L\n",
    "Yq4d_log = 20*np.log10(Yq4d_norm)\n",
    "\n",
    "fig = plt.figure(figsize=(10,4))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(n[0:L//2+1]*fs/L, Yq4_log[0:L//2+1])\n",
    "axes.plot(n[0:L//2+1]*fs/L, Yq4d_log[0:L//2+1])\n",
    "axes.set_title('Comparing Power Spectra with and without Dithering')\n",
    "axes.set_xlabel('Frequency (Hz)')\n",
    "axes.set_ylabel('Power Spectra (dB)')\n",
    "axes.grid(True, 'Major')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that dithering has broken-up periodic elements in the quantised signal, causing the spurious tones to be suppressed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Conclusion\n",
    "This notebook explored the quantisation process and has introduce dithering.\n",
    "\n",
    "The next notebook investigates ADCs and DACs."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "[⬅️ Previous Notebook](01_sampling.ipynb) || [Next Notebook 🚀](03_adcs_and_dacs.ipynb)\n",
    "\n",
    "Copyright © 2023 Strathclyde Academic Media\n",
    "\n",
    "---\n",
    "---"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
